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Abstract. A simple one-dimensional cellular automaton model with threshold 
dynamics is introduced. It is loaded at a uniform rate and unloaded by abrupt 
relaxations. The cumulative distribution of the size of the relaxations is analytically 
computed and behaves as a power law with an exponent equal to —1. This coincides 
with the phenomenological Gutenberg-Richter behavior observed in Seismology for the 
cumulative statistics of earthquakes at the regional or global scale. The key point of 
the model is the zero-load state of the system after the occurrence of any relaxation, 
no matter what its size. This leads to an equipartition of probability between all 
possible load configurations in the system during the successive loading cycles. Each 
cycle ends with the occurrence of the greatest -or characteristic- relaxation in the 
system. The duration of the cycles in the model is statistically distributed with a 
coefficient of variation ranging from 0.5 to 1. The predictability of the characteristic 
relaxations is evaluated by means of error diagrams. This model illustrates the value of 
taking into account the refractory periods to obtain a considerable gain in the quality 
of the predictions. 
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1. Introduction 

The cellular automaton model presented here will be denoted, for brevity, as the C- 
Model (CM). It has analogies and differences with the so called Minimalist Model 
(MM) [2] that will be emphasized in the next Sections. These types of cellular 
automaton models appeared in the context of the Self Organised Criticality (SOC) 
paradigm, introduced by Bak et al. [3^ .4], which meant a very important new 
conceptual perspective to try to understand complexity in nature in general and in 
natural hazards, such as forest-fires, landslides and earthquakes, in particular. The 
flagship of SOC is the sand-pile model. It is a conservative cellular automaton where 
the size-frequency distribution of avalanches exhibits a power-law behaviour where the 
exponent corresponding to the non-cumulative distribution is around —1. 

We will pay particular interest to studying two properties of the CM, namely 
the probability of occurrence of a relaxation of size k for a system of size A^, and 
the probability distribution for the time interval between successive characteristic 
relaxations, i.e. for the length of the loading cycles. This distribution will be called 
PnIu), where n is the discrete time elapsed since the occurrence of the last characteristic 
relaxation. The mean and the variance of P/v(n) will receive special consideration for 
the reasons that will be explained later. 

For a better perspective of the results of this model we will refer to two important 
concepts coming from seismology. The best established law in regional seismicity is the 
Gutenberg-Richter relation between the magnitude of an earthquake and its frequency. 
This law is of the power-law type when magnitudes are expressed in terms of rupture 
areas 



to or greater than S, and b is the so-called 6- value, which is a universal phenomenological 
parameter with a value close to unity [5] . Note that equation ([T]) represents a cumulative 
distribution. The corresponding non-cumulative distribution would also be a fractal law 
with an exponent equal to 



The interpretation of 6 as a sort of universal critical exponent in nature exactly equal to 
1 has spurred theoreticians for years and many ideas, typically based on specific models 
and/or mechanisms, have been offered to explain both the power law behavior of ([1]) 
and the value of the b parameter P, [TJ El [9l [ini [H] • Therefore, the Sand Pile model 
is not suitable for describing the observed Gutenberg-Richter Relation. But other SOC 
cellular automaton models for earthquakes, such as those introduced by Olami, Feder 
and Christensen ^ provided, within a certain limit of conservation, a correct order of 
magnitude for the b parameter. For a review of SOC and earthquakes, see [121 US]- 

Besides, it is important to bear in mind that the above mentioned Gutenberg- 
Richter law is a property of regional seismicity, appearing when one averages seismicity 





(2) 
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over big enough areas and long enough time intervals (e.g. [S]). For individual faults, 
the type of size-frequency relationship is different from the G-R law and consists of an 
approximate power-law distribution of small events (small compared with the maximum 
earthquake size that a fault can support, given its area), which occur between roughly 
quasi-periodic earthquakes of much larger size that rupture the entire fault. These 
large quasi-periodic earthquakes are termed "characteristic" [15], and the resulting 
size-frequency relationship is the Characteristic Earthquake distribution. It must be 
mentioned that the very concept of characteristic earthquake is under debate [T6t fT7]. 

Although the terminology that we use in this paper is often reminiscent of the 
seismological jargon, the CM could be applied to any entity where energy enters at 
a constant smooth rate and exits in the form of sudden relaxations with a very low 
efficiency because of dissipation. 

The structure of the paper is as follows. In Section [21 the rules of the model are 
given together with a table and several figures of numerical results. Section [3] is devoted 
to some analytical deductions. In particular, using several properties of Markov chains, 
we compute the stationary probabilities of the configurations in the CM, the probability 
of occurrence of relaxations of any size, and some other properties of the loading cycle. 
Section mis devoted to the evaluation of the predictabihty of the large relaxations in this 
model. This is accomplished by using the so called error diagrams. A brief discussion 
and the conclusions are gathered in Section [5l In a final Appendix we explain some 
detailed calculations for the case = 3. 

2. The C-Model: rules and some numerical results 

The CM is a simple one-dimensional cellular automaton with threshold dynamics. 
Consider a one- dimensional array of length A^. The ordered positions in the array 
will be labeled by an integer index i varying from 1 to A^. This system is loaded 
by receiving individual particles in the various positions of the array, and unloaded by 
emitting groups of particles through the first position i = 1, which are called relaxations 
or events. These two functions proceed using the following five rules: 

(i) The incoming particles arrive at the system at a constant rate. Thus the time 
interval between successive particles will be the basic time unit in the evolution of 
the system. In comparison with this time interval, the time taken by the relaxations 
is negligible. 

(ii) All the A^ positions have the same probability of receiving a new particle. When a 
position receives a particle we say that it is occupied. 

(iii) The reception of a particle saturates a position, i.e. if a new particle comes to a 
position already occupied, that particle is dissipated. 

(iv) The position i = 1 is special. When a particle goes to the first position a relaxation 
occurs. Then if all successive positions from i = 1 to i = k are occupied and 
the position A; + 1 is empty, the effect of the relaxation is to unload all the levels 
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Figure 1. (a) Size-frequency plot for tlie relaxations in the CM; the three cases 
N — 10, 100, and 1000 are superimposed, (b) cumulative magnitude-frequency plot. 
In both graphs a line of slope —2 (top) and —1 (bottom) has been added for visual 
comparison. 



from i = 1 to i = k. Hence the size of the event is k. The maximum relaxation 
corresponds to k = N, which is called the characteristic relaxation of the system. 

(v) Besides, those positions that were occupied by particles but were not affected by 
the relaxation lose their particles so that the array is left completely empty. 

Thus in the CM, after any small relaxation, the system always restarts from the 
empty state, keeping no memory. The simplicity of these five rules makes it very 
appropriate for numerical simulations. The formalism and results of the Markov Chains 
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Table 1. Probability of occurrence of a relaxation of magnitude k for a system size, 
N. 



Size 


N=2 


N=3 


N=4 


N=10 


N=100 


k=l 


0.50006 


0.50003 


0.50000 


0.50009 


0.49948 


k=2 


0.49994 


0.16666 


0.16652 


0.16647 


0.16667 


k=3 


- 


0.33331 


0.08333 


0.08348 


0.08359 


k=4 


- 


- 


0.25015 


0.04994 


0.05012 


k=5 








0.03339 


0.03358 


k=6 








0.02370 


0.02370 


k=7 








0.01785 


0.01768 


k=8 








0.01396 


0.01413 


k=9 








0.01108 


0.01118 


k=10 








0.10002 


0.00906 


k=99 










0.00009 


k=100 










0.00990 



can also be applied to this model. In contrast with the CM, in the MM [Ij, rule 5 did 
not exist. The fifth rule of this model could be representative of what occurs with the 
precipitation of a cloud. In the process of growth of the CCN (cloud condensation nuclei) 
by condensation, only those CCN that have been activated become cloud droplets and 
take part in the posterior processes of collision and coalescence, and therefore only they 
eventually contribute to the rain. The so-called haze droplets that were not activated 
evaporate. 

Using these five rules and carrying out numerical simulations, we have obtained 
the results contained in Table [T] and in figures [TJHl and figures [6] and lAll later. From 
Table [Hand Fig. [TK we see that in the CM the value of pk{k < N) is an iV-independent 
constant. This sort of scale invariance, which also appeared in the MM, can be justified 
by using the same arguments written in Ref. [1]. We also note that the non-cumulative 
size-frequency plot for the events in this model is of the Characteristic Earthquake type, 
as mentioned in the Introduction. In Fig. [T)d, it is shown that, independently of N, the 
cumulative size-frequency plot for the relaxations in this model is a perfect power law 
with an exponent equal to —1. It is worth saying at this point that, as far as the 
authors know, there is no other model with this beautiful property. In the next Section, 
a probabilistic argument will be provided for this notable result. 

Note that the non-cumulative distribution in Fig. [Tti clearly shows the special 
character of the maximum magnitude relaxation, while the cumulative distribution hides 
this preponderance. 

In Fig. [2], we show the probability distribution for the length of the cycles, P/v(n), 
for = 3 (squares) and = 10 (circles). In both cases the distribution shows an 
initial refractory period (i.e., an interval of null probability of relaxation occurrence) of 
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Figure 2. Probability distribution for the length of the cycles, PN{n), for TV = 3 
(squares, left axis), and = 10 (circles, right axis). In the iV = 3 case the simulation 
result is compared with the analytical result (continuous line). 



length equal to A^ — 1, and a slow asymptotic decrease which will be reflected in the high 
variances of this model. In case Psin), the results obtained by simulations are compared 
with the analytical result (continuous line) obtained in the Appendix. 

3. Stationary Probabilities in the CM, Spectrum of Earthquakes and other 
Analytic Properties 

In the CM, for a system of size A^, there are N groups of configurations distinguished 
by their index of load occupation 6, 

9 = 0,1,2,...,N -1. (3) 

In each of these groups there are 



N-1 

e 



(4) 



different stationary configurations denoted by 9j, where the integer index j runs in the 
range 

As in any Markov chain of this type, if the Markov matrix is denoted by M, then the 
stationary probabilities for each configuration are the components of the eigenvector of 
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M-^ corresponding to the eigenvalue +1. These components will be denoted by i^{9j). 
Their j-independent values are 

^(^.) = ^ / X • (6) 



A^- 1 



Relations ([6]) are identified in a straightforward way because of the regularity of the 
Markov matrix M in this model (as an illustration, the case = 3 is worked out in 
detail in the Appendix). Thus, all the configurations that share the same 9 have the 
same probability, and the sum of the probabilities of all of them is 1/A^. In other words, 
any load 9 in this system has identical probability 1/A^, and the probability among all 
the possible configurations corresponding to a given 9 is also equally distributed. 

Using these results, let us compute first the fraction of particles that, on their 
way to a position in the array, find this position occupied and are therefore lost. This 
proportion will be called the refiectivity p of the system: 

9 




-Jp 1^ ^- m 9 -o(^-n}- 

0=0 

Thus, for large A^ nearly 50% of the particles are lost by reflection. This result is 
compared with the simulations in figure [31 We see there that the refiectivity (open 
and filled squares) increases with the system size, abruptly at the beginning and more 
slowly for bigger system sizes, reaching an asymptotic value of ~ 0.5 for big systems 
(A^ > 200). The smallest refiectivity is p ~ 0.34 for A^ = 3. 

A second analytical conclusion is the following. From Eq.(l6]), we know that in the 
case of maximum occupation, i.e., when there are A^ — 1 particles in the system, its 
probability is 

vr(iV-l) = ^, (8) 

therefore the probability that at any arbitrary time a characteristic event occurs is 
1/A^^ and, in consequence, the mean value of the time elapsed between two consecutive 
characteristic relaxations in this model is 

< P^in) >= N\ (9) 

This conclusion coincides with the numerical result presented in figure H^, better seen 
in the inset where log-log scales have been used to appreciate the +2-slope trend of 
the mean. Note also how, for small system sizes, the mean is bigger than the standard 
deviation but that both of them tend to a common value as the system size increases. 
This can be better appreciated in figure |Dd where the aperiodicity, i.e. the ratio of the 
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System size, N 



Figure 3. Reflectivity p, efficiency e, and fraction of lost particles during relaxations 
A as a function of N. For details, see Section [3l 

standard deviation to the mean, is plotted. The minimum aperiodicity is 0.57 (A^ = 3) 
and grows asymptotically to unity. 

We also know that the mean time between two arbitrary consecutive relaxations in 
the model is N , and therefore we deduce that in this model, between two consecutive 
characteristic events, there are, on average, N — 1 non-characteristic relaxations. 

As mentioned above, a cycle is the time between two successive characteristic events. 
And each cycle is formed by a succession of several initial sub-cycles, marked by the 
occurrence of non characteristic events, and a final sub-cycle marked by the occurrence 
of the characteristic relaxation that ends the cycle. The calculation of the length of the 
final sub-cycle is straightforward. This time length is formed by the addition of two 
independent stages. The first stage, which starts from the state 6' = and ends when 
the system arrives to the state 9 = N — 1, will be called the loading stage. The second, 
in which the system is completely loaded, ends when a particle hits the first site; this 
will be called the hitting stage. The length of the loading stage is the mean time taken 
to fill a box of A^ — 1 cells by random assignments of the particles taking into account 
the reflection of particles by occupied cells [18]. The length of the hitting stage is the 
mean time of a geometric distribution with a hitting probability equal to 1/A^. Thus, 
we have 



< n >fsc=< n>i + <n >h= N{S - 1) + N, 



(10) 
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System Size, N 



Figure 4. (a) Mean and standard deviation of PN{n) as a function of N\ the inset 
shows the trend of the mean in log-log scales. MC stands for MonteCarlo results, (b) 
Aperiodicity versus N . 



where 5* is a function of N defined as 

N 



k=l 



For > 10, S" is already well approximated by 

lim S = C + \tiN + — + o(- , 

N^oo 2N \ m 



\im S = C + \tiN + \ + o(^] , (12) 
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where C = 0.577215 ... is Euler's constant. Therefore, the asymptotic behavior of the 
duration of the final sub-cycle is of the form 

lim < n >fsc= NlnN. (13) 

The variance of the time length of the final sub-cycle can also be analytically 
calculated. 

Using numerical simulations, we have also seen in Section [2] that the CM provides 
a power-law relation, with 6=1, for the cumulative size-frequency distribution of 
the relaxations. This power-law relation can be understood by a simple probabilistic 
argument. In the CM, after any relaxation, the grid contains no particles, i.e., all the 
positions of the system are empty. Therefore, to have in the next event a relaxation of 
size > k {1 < k < N), it is necessary that the process of occupation of the lowest k 
levels be in such a way that the first level {i = 1) is the last one to be occupied. The 
question then is: What is the probability that k identical sites be occupied with the 
condition that one of them be the last? As all of the sites are equivalent, the answer is 

F 

Thus, for 1 < A; < AT, 

P{> k) = l (15) 

independently of the size of the system. 
Now, due to the fact that 

p{k) =p{>k) -p{>k + l), (16) 

which can be applied for 1<A;<A^ — l,we find that 

^^^^ ^k~ ITl ^ k{k + l) ^ F(l + 1/A;)' 

k<N~l. (17) 



For the limit case k = N, we apply expression (1151) and obtain 

p(k = N) =p{k>N) = ^. (18) 
The verification of the adequate normalization of p{k) is carried out as follows: 

N-l N-l 



/I 1 \ 1 

I — 1 L — 1 V / 



k=l k=l 

Therefore, in the CM the cumulative spectrum of relaxation given in (11 5p constitutes a 
power law with an exponent —1. This coincides with the phenomeno logical Gutenberg 
- Richter exponent observed in seismology for the cumulative statistics at the regional 
or global scale. And according to this, in the CM the non-cumulative form of the size- 
frequency relation, as given in (fTTll . is asymptotically a power law with an exponent 
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(20) 



equal to —2. These conclusions exactly reproduce the plots in figured! Relation (ITBI) is 
the discrete derivation of f|T5|) : in the continuum it would correspond to 
d I _ I 

With respect to the probability partition that relates a system of size with the 
next one of size + 1, that is, the form in which relation 

P;v(iV) =P7V+l(iV)+P7V+l(iV + l) (21) 

is implemented in this model, we have the following simple form 

iV ~ iV(iVTI) ^ N +1' 

Here, it is interesting to note that the counterpart of equation ( 122|) in the MM can only 
be obtained numerically. 

Let us now compute the efficiency of the system e, i.e., the proportion of particles 
that contribute to the relaxations of the system. We know that, on average, in time 
units a relaxation is produced. Thus, 

[N-l 



1 



(23) 



And using the function S'(A^) as defined in (fTTj) . the efficiency can be expressed as 

e = -. (24) 

This function amounts to e = 0.61 for A^ = 3 and monotonously decays to as A^ grows 
(figure [3], open and filled stars). Thus, we realise that this model is by far less generous 
than a standard slot machine. For example, present laws in Nevada impose a minimum 
of e = 0.75 for the payout percentage of slot machines in the casinos of this state. For 
large A^, 

lim e = (25) 

In this stochastic model, the history of one particle can be only one of the three following 
things: (i) refiected by the system, (ii) emitted in a relaxation, or (iii) lost on the occasion 
of a relaxation in the system. Calling A to the fraction of these lost particles, we have 

A = 1 - p - e. (26) 

Due to the fact that e ^ for large A^, in this limit both p and A share one half of the 
probability. This is illustrated in figure [3] where we see how A (open and filled circles) 
grows from A ~ 0.05 for A^ = 3 to A = 0.5 for big system sizes. 
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4. Forecasting the Large- Avalanche Occurrence. 

A hint of the predictability of the large relaxations in this type of model is given by 
the aperiodicity of their time series. The aperiodicity is a quantitative measure of the 
lack of regularity of a time series. If fi is the average time between two consecutive 
characteristic relaxations (i.e., the mean duration of the cycle), and a is the standard 
deviation of the duration around the mean, then a = a/fi. The aperiodicity is otherwise 
known as the coefficient of variation. A value of a = means that the system is periodic 
and all the cycles have exactly the same duration; the range < a < 1 defines quasi- 
periodic time series; and the case a = 1 can correspond to a purely random (Poissonian) 
time series, but not necessarily so. Time series with a > 1 are said to be clustered. As 
the a of this model takes values between 0.5 and 1, the occurrence of the large events is 
a quasi-periodic phenomenon. A robust way to assess the predictability of a time series 
is by trying to forecast its events by declaring alarms at particular times (Figure [5]). 
The aim is to declare alarms before all the events in order not to miss any event, but 
to declare them just before the events in order to minimize the total alarm time. Many 
strategies can be devised to declare the alarms but there is a reference strategy to which 
all others can be compared [2l|T9l[20]. This strategy consists of setting the alarm a fixed 
time interval after each event (waiting time, t) and maintaining it until the occurrence 
of the event. If the following event in the time series occurs before the alarm is raised, 
it is counted as a prediction error; if the following event in the time series occurs after 
the alarm is raised, it is counted as a prediction success and the alarm is then cancelled. 

The fraction of errors /e (number of missed events divided by the total number of 
events) and the fraction of alarm time fa (total alarm time divided by the total duration 
of the time series) can be computed as a function of the above mentioned waiting time 
t, and the purpose is to find the optimum waiting time. This optimum waiting time 
depends on the relative importance that failing to predict an event has compared to 
keeping the alarm on. An objective function, called loss function, L, can be defined 
that incorporates this trade-off in each particular case. Here we will use the simplest of 
them, L = fe + fa, where failure to predict and a long alarm time are equally penalized. 
Thus, our aim is to find the waiting time t = t* that minimizes L{t). This minimum 
value is denoted by L* = L{t*). And the best way to graphically display this is by means 
of an error diagram, where the fraction of errors /e runs along the horizontal axis and 
the fraction of alarm time fa runs along the vertical axis (Fig. [5]). Error diagrams were 
introduced in earthquake forecasting by Molchan [21] who contributed with rigorous 
mathematical analysis to the optimization of the earthquake prediction strategies. 

Figures y Wp show the results of applying the reference strategy to a system of 
A = 10 and A = 100 . For A = 10 : t* = 29, = 0.09, /„ = 0.72 and L* = 0.81; 
for A = 100 : t* = 588, fe = 0.02, fa = 0.94 and L* = 0.96. This increase in L* is 
what one would reasonably expect as the aperodicity of these systems is a = 0.78 for 
A = 10 and a = 0.96 for A = 100. The question now is, can we perform a better 
job in forecasting the occurrence of the large relaxations in this model? The answer is 



13 

ESS S 

12 3 4 

waiting 
time 




Figure 5. Reference strategy for the assessment of the predictabihty of a time series 
and its representation on an error diagram. The events that are to be predicted (large 
relaxations) are the vertical red bars numbered correlatively. An alarm is set a fixed 
time interval after each event (waiting time) and the prediction is labeled error (E) or 
success (S) depending on whether the alarm was off or on when the event occurred. 
The fraction of errors is the number of events not predicted (one in the example) 
divided by the total number of events (four events), i.e., fe ~ 0.25; and the fraction 
of alarm time is the total alarm time (blue sections of the time line, 71 time units) 
divided by the total duration of the time series(178 time units), i.e., fa = 0.4 in the 
example shown in the figure. 

definitely yes. For this goal we will exploit the conspicuous refractory periods that take 
place after the occurrence of any relaxation in the model. 

The reference strategy consists in connecting the alarm at t time units after the 
occurrence of each characteristic relaxation in the system, and the best t is that which 
minimizes the loss function L. As the alarm is not disconnected until the occurrence of 
the next characteristic relaxation, in this strategy there are no false alarms. The new 
strategy consists in connecting the alarm at {N — l)+t' time units after the occurrence 
of any relaxation, big or small, in the system. And the best t' is that which minimizes 
L. As in this strategy the alarm is disconnected when any relaxation takes place, we will 
have false alarms that will contribute to the total fa- The results of the new strategy 
are plotted in figures [6t and [6li. The best values of the parameters for = 10 are 
t'* = 7, /e = 0.05, fa = 0.18 and L* = 0.23; and for A^ = 100 are t'* = 214, /e = 0.01, 
fa = 0.04 and L* = 0.05. It is clear that the new strategy has been a success. 

5. Discussion and Conclusions 

The CM presented here is a new model of the type that started with the MM |Tl |2]. In 
these models, the only parameter to be fixed is A^, the total number of sites, or cells, 
susceptible to occupation. The two models are materialized on discrete one- dimensional 
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Figure 6. Error Diagrams (a)(c) for = 10 and (b)(d) for N — 100. The results of 
the Reference Strategy are shown at the left. The right figures show the improvement 
in the predictability provoked by considering the refractory periods. For details see 
the text. 



arrays, which permits their study by means of Markov chains or straightforward 
Montecarlo simulations. 

The CM is by no means conservative: the load is dissipated in several forms. The 
array that is loaded and unloaded in successive cycles dissipates load by reflection of the 
already occupied sites, and when relaxations are produced. As a result of its rules, this 
system maintains a statistical equilibrium such that all possible values of the load in the 
system have the same probability. Besides, all the configurations with the same load are 
also equally probable. The non-cumulative size-frequency relation for all the relaxations 
in this model is of the Characteristic-Earthquake type, with an excess of system-wide 
events. After an initial transient, the probability of the small relaxations is a power 
law with an exponent equal to —2 (see figure [TK). In consequence, the cumulative size- 
frequency relation for all the events in this model is a neat power-law with exponent 
— 1. This is notable because it coincides with the form of the Gutenberg- Richter law of 
the global seismicity. 

The lack of memory in the CM with respect to the MM makes this model simpler. 
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A clear manifestation of this simplicity can be observed in the computation of the 
probability distribution of the length of the cycles when the Jordan form of the respective 
M' matrices are computed in the CM (see the Appendix for the definition of matrix M' 
and other algebraic details). In the CM the Jordan form of M' is a diagonal matrix and 
therefore PAr(n) is just a combination of exponentials (with n in the exponents). In the 
MM, however, the appearance of off-diagonal terms in M' led PN{n) to be formed by 
terms where polynomial terms in n are multiplied by exponential terms in n |22j . 

In spite of the fact that in the CM the period of return of the big relaxations has a 
aperiodicity between 0.5, for small N, and tends to 1, for large N, we have shown that 
these big relaxations can be accurately predicted by exploiting the refractory periods 
that take place after the occurrence of any relaxations in the system. In fact, it has 
been shown that, the bigger N is, the better can be the result of the forecasting. 

There are plenty of natural systems where just after the occurrence of a 
characteristic event, the probability of occurrence of a new one is null or very low. 
Examples of this type of inhibition are the process of firing of neurons, the refractory 
periods that occur between the breakouts of contagious diseases, etc. This model has 
illustrated the importance of taking into account these periods for making accurate 
predictions. 

5.1. Note Added. 

An anonimous referee has suggested us to make a comparison between the model 
studied in this work, the CM (and especially the MM) and those models recently 
used to describe the dynamics of microtubules. Microtubules are linear polymers that 
serve as structural components within cells and are involved in many cellular processes 
including vesicular transport. In its simplest version, called the minimalist model [23] 
microtubules are represented by a chain formed by two types of monomers: GTP 
(for guanosine triphosphate) and GDP (for guanosine diphosphate). The state of a 
microtubule evolves due to the three following three processes: 

(i) Attachment. A microtubule grows by attachment of a GTP monomer at its tip, 
with a rate A if in the tip there is a monomer GTP, and with a rate p\ if in the tip 
there is a GDP monomer 

(ii) Conversion. Each GTP monomer of the chain can be independently converted by 
hydrolysis into a GDP monomer. This occurs at a rate unity. 

(iii) Detachment. A microtubule shrinks due to the detachment of a GDP monomer 
that is situated at the tip of the microtubule. This process takes place at a rate . 

In the phase space subtended by the set of parameters (A, /x, p) there is a rich phe- 
nomenology for the microtubule growth dynamics such as phase boundaries separating 
regions where the microtubule grows, on average, at a certain rate from other regions 
where the average length of the microtubule is finite, etc. And as a glaring event in the 
limit of large /i it is worth mentioning that a global catastrophe can occur when a newly 
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attached GTP at the tip converts to a GDP and the rest of the microtubule consists 
only of GDP at that moment, the microtubule instantaneously shrinks to zero length. 

And the question now is which are the differences and similarities between these 
idealized microtubule growth models (MGM) and our CM and MM? Both type of mod- 
els are stochastic and in one spatial dimension but while the MGM describe a system 
where its length can grow or shrink, i.e. fluctuates in time, in both the CM and the 
MM the length of the array where the particles are positioned, A^, is fixed. Obviously 
this difference is crucial. What fluctuates in our models is the number of particles, k, 
emitted in a relaxation, which varies between 1 and A^. 

Having stated the main difference, the question would be: could we modify our 
models to resemble the MGM? Maintaining all the other rules of these models, this task 
would require addding new basic ingredients. Assuming that the flrst site corresponds to 
the tip of the polymer, a new parameter a with dimensions of a rate would be responsible 
of adding new free sites at the tip. Second, after the occurrence of a relaxation of k 
particles, the k sites of the array neighbour of the tip would disappear and the k + 1 
position of the array, which is free, would turn into the new flrst site responsible of 
controlling the new relaxation. Thus, in a short dictionary between the two models a 
GTP monomer would be represented by a free site, a GDP monomer by an occupied 
site, the reception of a particle shifts the site from GTP to GDP, a relaxation shrinks 
the length, and characteristic relaxation would be the equivalent to the catastrophe of 
the microtubule. The new MG-like model sketched here would be more economical in 
parameters than those studied in [23]. Its properties will be studied in a next future. 

Appendix A. Explicit calculations for A^ = 3 

In this Appendix the algebraic calculations corresponding to the case A^ = 3 are worked 
out in detail. These results will likely be useful to the reader for a better understanding 
of the content of Section [3l 

A scheme of the possible states for this system is plotted in Fig. IA1[ Note that 
the relation between this notation and the more general notation expressed in equations 
©, Q and dSD is a = = 0), b = {9 = 1, j = 1), c = {9 = 1, j = 2),d = {9 = 2). 

Using the rules of the model, the probabilities of the non-null one-step transitions 
are: M,^, = 1/3, M,^, = 1/3, M,^, = 1/3, M,^, = 1/3, M^^b = 1/3, M,^d = 1/3, 
M,^, = 1/3, M,^e = 1/3, M,^, = 1/3, M^^a = 1/3, M.^^ = 2/3. Thus, for A^ = 3, 
the Markov matrix is 



/ 1/3 


1/3 


1/3 


\ 


1/3 


1/3 





1/3 


1/3 





1/3 


1/3 


V 1/3 








2/3 J 



(A.l) 



The stationary probabilities of residence in this Markov system are the normalized 
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b c 
Configuration 



Figure Al. Configurations and residence probabilities for the case iV = 3. Only four 
states, labeled a, b, c, and d, are possible (inset). The probability of states a and d is 
1/3 and the probability of states b and c is 1/6. 



components of the eigenvector of corresponding to the eigenvalue A 
the transpose of M. These components are 

/ 1 \ 



lA 



1/2 
1/2 

V 1 / 



1, being 



(A.2) 



which agrees with that written in (j4]) and ([6]) for any A^. In the calculation of the 
probability distribution of the length of the cycles we will follow the general method 
used in [T], that is 

1 



PNin) 



N 



[M 



nn-li 



0,9 = N-1] 



(A.3) 



This equations means that the probability of having a cycle of length n is given by 1/A^ 
times element = 0, ^ = — 1) of the matrix [M']"~^. This is justified as follows. 
On the left, we have, by definition, the probability that the cycle lasts n steps. On the 
right, it is written the probability that, in n — 1 steps, the system has moved from the 
empty state to the state with a load A^ — 1, with the restriction that this is the first visit 
to this maximally loaded state. This restriction is easily implemented by making the 
matrix element (6' = A^ — 1, ^ = 0) in M equal to 0. That "pruning" converts matrix M 
into the new matrix M'. Once the system has done this transit in n — 1 steps, with the 
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mentioned restriction, the probability of passing in one step from the state ^ = — 1 
to the state 6' = is This is the reason to include this factor in flA.3l) . 

In the case N = 2 the MM p] and the CM are actually the same model. The P{n) 
distribution is 

(A.4) 



PN=2{n) = 
In the case = 3, M' is 

/ 1 1 1 \ 



M' 



(A.5) 



Q 



v 



rn-l 



1 

1 

n-1 

X 



v 
















4^2 



Inserting the corresponding formulae into (lA.Sp . we obtain 

PN=3in) = 



2(1 - v^) 



2 + V2 
2(1 + v^) 



(l + v^)"Ln>2. 



(A.6) 
(A.7) 



1 110 1 
3 10 11 
^ 2 / 

In order to obtain a generic power of this matrix, we will perform its Jordan 
decomposition 

M' = QJQ-^ 

and hence 

M'"^ = qr^Q \ 

where the factors on the right of (]A.7p are 

/ -2 -^2 V2 \ 
-1-11 1 
-1 1 1 

y 

■r " ' 

( 1 


(1 + 72)"-^ j 

( -1/2 1/2 \ 
1 

-iTl V4 1/4 
V 1/4 1/4 ---^ 



(A. 



(A.9) 
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This function has been plotted in figure [2] together with its numerical counterpart. It 
is null for n = 2, as it should, because for = 3 the minimum length of a cycle is 3. 
It is apparent that PN=3{n) is formed by a combination of three geometric-like terms. 
This is due to the fact that matrix J has no off-diagonal term. Inside the brackets 
the first term is negative, the second is oscillating, and the third -and dominant- is 
positive. Therefore, in the limit of large n this function is composed only by the third 
positive geometric decaying term and, in consequence, the asymptotic hazard rate of 
this model is a constant. The mean fi and standard deviation a of PN=3{'n) are 9 and 
271/2 respectively, with an aperiodicity a = a / ji = 0.58. 
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